**** COX ET AL. (2023) ADAPTATION OF RAMEY/ZUBAIRY CODE. 
**** This code and the underlying data file (rzdat.xlsx) comes directly 
**** from the replication package of Ramey Zubairy (2018). 
**** We adapt it slightly to estimate their IRFs using sub-combonents of G 
**** rather than G as a whole. 


**** JORDAGK.DO

*** May 31, 2016
*** This program estimates IRFs and both 2-step and 1-step multipliers using the Jorda method

*** Uses Gordon-Krenn (GK) transformation (i.e. variables are divided by potential)

***
*** Requires:
***     rzdat.xlsx, updated April 7, 2016
********************************************************************************

#delimit;

drop _all;
clear all;

set more 1;
set matsize 800;

capture log close;
log using jordagkirfs_results.log, replace;


/*******************************************************************************
  SET PARAMETERS THAT GOVERN SPECIFICATION
*******************************************************************************/

local sample = 1;  /*1 = full sample, 2 = post-WWII */

local omit nomit;  /*either nomit(don't omit subsample) or wwii (omit WWII) */

local state slack;  /* slack or zlb or recession or ag*/

local shock newsy; /* shock identification: either newsy or bp */

local p = 4; /*number of lags of control variables*/

local trends = 0; /*0 = no trends, 1 = trends */

local tax = 0; /*0 = exclude taxes, 1 = include taxes */


*******************************************************************************;
** RAW DATA IMPORTATION AND DATA SETUP;
*******************************************************************************;

import excel ../data/rzdat.xlsx, sheet("rzdat") firstrow;

*insheet using rzdatnew.csv;

drop if quarter<1889;

gen qdate = q(1889q1) + _n-1;
tsset qdate, q;

// Merge in contracts data (COX ET AL ADAPTATION) 
merge 1:1 quarter using "../data/contracts_for_ramey_merge.dta";


/* World War II rationing sample.  Start in 1941q3 because of constraints on
     auto production.  Most rationing ended in Aug/Sept 1945, a few items in
	 Nov/Dec 1945 */

gen wwii = quarter>=1941.5 & quarter<1946;

gen nomit = 0;  /* indicator for no omit */


*** DEFINE QUARTIC TREND;

gen t = _n;
gen t2 = t^2;
gen t3 = t^3;
gen t4 = t^4;

*** DEFINE STATE VARIABLE;

gen slack = unemp >= 6.5;  /* unemployment state with fixed threshold */
gen slack8 = unemp>=8;
gen slackhp = unemp>=hpunemp_split;  /* unemployment state with hp threshold */

gen zlb = zlb_dummy;  /*  zlb state */
gen zlb5 = tbill<=0.5; /*tbill rate less than 0.5*/
gen both =unemp >= 6.5 & zlb_dummy>=1; /* consider slack and ZLB state together*/

*** NORMALIZATION;

/* choice of potential GDP for normalization:

   rgdp_potCBO (cubic trend early, CBO late) or rgdp_pott6 (6th degree for full),
   both fitted excluding Great Depression, WWII:  quarter>=1930 & quarter<1947*/

local ynorm rgdp_pott6; /* rgdp_pott6 or rgdp_potcbo */

sort qdate; 

* BASIC VARIABLES;
gen newsy = news/(L.`ynorm'*L.pgdp);
gen rgov = ngov/pgdp;

/* New Variables */ 
gen rcontracts = contract_proxy/pgdp; 
gen rwages     = wages/pgdp;
gen rnon_wage  = non_wage/pgdp; 

gen rtax = nfedcurrreceipts_nipa/pgdp;
gen taxy = nfedcurrreceipts_nipa/ngdp;
gen debty = pubfeddebt_treas/L.ngdp;
gen lpgdp = ln(pgdp);
gen ly = ln(rgdp);

gen infl = 400*D.lpgdp;

* normalize variables and shorten names;

gen y = rgdp/`ynorm';
gen g = rgov/`ynorm';
gen c = rcontracts/`ynorm'; 
gen w = rwages/`ynorm';
gen n = rnon_wage/`ynorm';
 
gen bp = g; /* Blanchard-Perotti shock is just orthogonalized current g */

*** AG DEFINITION OF STATE:  ag = 1 is extreme recession, ag = 0 is extreme expansion;

gen z = 100*(F3.ly - L4.ly)/7;  /* AG definition of state */

*The mean of z is approx. 0.8 and std is 0.5.  AG specically use those numbers so we do too;

gen znorm = (z - 0.8)/0.5;

gen fznorm = exp(-1.5*znorm)/(1 + exp(-1.5*znorm));

gen ag = fznorm;

*******************************************************************************;
** CUMULATIVE VARIABLES;
*******************************************************************************;

gen cumuly = 0;
gen cumulg = 0;
gen cumulc = 0; 
gen cumulw = 0; 
gen cumuln = 0; 
 
forvalues i = 0/20 {;

   gen f`i'cumuly = F`i'.y + cumuly;
   gen f`i'cumulg = F`i'.g + cumulg;
   gen f`i'cumulc = F`i'.c + cumulc; 
   gen f`i'cumulw = F`i'.w + cumulw; 
   gen f`i'cumuln = F`i'.n + cumuln; 
   
   gen recf`i'cumulg = f`i'cumulg*L.`state';
   gen expf`i'cumulg = f`i'cumulg*(1-L.`state');
   
   gen recf`i'cumulc = f`i'cumulc*L.`state';
   gen expf`i'cumulc = f`i'cumulc*(1-L.`state');
   
   gen recf`i'cumulw = f`i'cumulw*L.`state';
   gen expf`i'cumulw = f`i'cumulw*(1-L.`state');
   
   gen recf`i'cumuln = f`i'cumuln*L.`state';
   gen expf`i'cumuln = f`i'cumuln*(1-L.`state');
   
   replace cumuly = f`i'cumuly;
   replace cumulg = f`i'cumulg;
   
};


*******************************************************************************;
**  INTERACTION OF SHOCKS WITH STATE;
*******************************************************************************;

 foreach var in newsy bp { ;
 
   gen rec0`var' = `var'*L.`state';
   gen exp0`var' = `var'*(1-L.`state');
 
 };

*******************************************************************************;
** CREATE LISTS;
*******************************************************************************;

   if `sample'==1 {;

         gen h = t - 1;  /* h is the horizon for the irfs */
         global trendlist t t2 t3 t4;
   };

    else {;
         drop if quarter<1947;
         gen h = t - 232 - 1;
         global trendlist t t2;
     };
	 
forvalues i = 1/`p' {; 

  foreach var in newsy y g c w n taxy debty infl{;

    gen rec`var'`i' = L`i'.`var'*L.`state';
    gen exp`var'`i' = L`i'.`var'*(1-L.`state');
 
  };
};

  if `trends'==0 {;
  
    if `tax'==0 {;
  
      global newsylinxlist L(1/`p').newsy L(1/`p').y L(1/`p').g; 
      global bplinxlist L(1/`p').y L(1/`p').g;
	  global newsynlxlist L.`state' recnewsy? expnewsy? recy? expy? recg? expg?;
      global bpnlxlist L.`state' recy? expy? recg? expg?;
	
	};
	
	else {;
	
	  global newsylinxlist L(1/`p').newsy L(1/`p').y L(1/`p').g L(1/`p').taxy; /*L(1/`p').infl;*/
      global bplinxlist L(1/`p').y L(1/`p').g L(1/`p').taxy; /* L(1/`p').infl;*/
	  global newsynlxlist L.`state' recnewsy? expnewsy? recy? expy? recg? expg? rectaxy? exptaxy?; /* expinfl? recinfl?;*/
      global bpnlxlist L.`state' recy? expy? recg? expg? rectaxy? exptaxy?; /* expinfl? recinfl?;*/
	
    };
  };
  
  else {;
  
    if `tax'==0 {;
	
      global newsylinxlist L(1/`p').newsy L(1/`p').y L(1/`p').g $trendlist;
      global bplinxlist L(1/`p').y L(1/`p').g $trendlist;
      global newsynlxlist L.`state' recnewsy? expnewsy? recy? expy? recg? expg? $trendlist;
      global bpnlxlist L.`state' recy? expy? recg? expg? $trendlist;
	
	};
	
	else {;
	
	global newsylinxlist L(1/`p').newsy L(1/`p').y L(1/`p').g L(1/`p').taxy $trendlist;
    global bplinxlist L(1/`p').y L(1/`p').g L(1/`p').taxy $trendlist;
	global newsynlxlist L.`state' recnewsy? expnewsy? recy? expy? recg? expg? rectaxy? exptaxy? $trendlist;
    global bpnlxlist L.`state' recy? expy? recg? expg? rectaxy? exptaxy? $trendlist;
	
    };
	
};


global newsylinshock newsy;
global newsynlshock rec0newsy exp0newsy;

global bplinshock bp;
global bpnlshock rec0bp exp0bp;


** INITIALIZE SUM OF EFFECTS TO 0 AND PARAMETERS SERIES TO MISSING;

gen sumliny = 0; gen sumling = 0; gen sumlinc = 0; gen sumlinw = 0; gen sumlinn = 0; 
gen sumexpy = 0; gen sumexpg = 0; gen sumexpc = 0; gen sumexpw = 0; gen sumexpn = 0; 
gen sumrecy = 0; gen sumrecg = 0; gen sumrecc = 0; gen sumrecw = 0; gen sumrecn = 0; 

foreach var in bylin byexp byrec bglin bgexp bgrec bclin bcexp bcrec bwlin bwexp bwrec bnlin bnexp bnrec up95bylin up95byexp up95byrec up95bglin up95bgexp up95bgrec up95bclin up95bcexp up95bcrec lo95bclin lo95bcexp lo95bcrec up95bwlin up95bwexp up95bwrec lo95bwlin lo95bwexp lo95bwrec up95bnlin up95bnexp up95bnrec lo95bnlin lo95bnexp lo95bnrec lo95bylin lo95byexp lo95byrec lo95bglin lo95bgexp lo95bgrec seylin seyexp seyrec seglin segexp segrec seclin secexp secrec sewlin sewexp sewrec senlin senexp senrec multlin multexp multrec {;
  
  quietly gen `var' = .;
  
}; 


*******************************************************************************;
** ESTIMATION OF IRFS
*******************************************************************************;

forvalues i = 0/0 {; /*Must treat horizon = 0 different in case the shock is BP*/

  ivreg2 F`i'.y $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bylinh`i' = _b[$`shock'linshock];
  
  gen seylinh`i' = _se[$`shock'linshock];
  
  ivreg2 F`i'.y exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen byexph`i' = _b[exp0`shock'];
  gen byrech`i' = _b[rec0`shock'];
  
  gen seyexph`i' = _se[exp0`shock'];
  gen seyrech`i' = _se[rec0`shock']; 
  
  
  if "`shock'" == "bp" {; 
  
  gen bglinh`i' = 1; 
  gen bclinh`i' = 1; 
  gen bwlinh`i' = 1; 
  gen bnlinh`i' = 1; 
  
  gen seglinh`i' = 0; 
  gen seclinh`i' = 0; 
  gen sewlinh`i' = 0; 
  gen senlinh`i' = 0; 
  
  gen bgexph`i' = 1; 
  gen bcexph`i' = 1; 
  gen bwexph`i' = 1; 
  gen bnexph`i' = 1; 
  
  gen bgrech`i' = 1;
  gen bcrech`i' = 1; 
  gen bwrech`i' = 1;
  gen bnrech`i' = 1; 
  
  gen segexph`i' = 0;
  gen secexph`i' = 0; 
  gen sewexph`i' = 0; 
  gen senexph`i' = 0; 
  
  gen segrech`i' = 0;
  gen secrech`i' = 0; 
  gen sewrech`i' = 0;
  gen senrech`i' = 0; 
  
};

else {;

  ivreg2 F`i'.g $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bglinh`i' = _b[$`shock'linshock];  
  gen seglinh`i' = _se[$`shock'linshock]; 

  ivreg2 F`i'.g exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen bgexph`i' = _b[exp0`shock']; 
  gen bgrech`i' = _b[rec0`shock'];
  
  gen segexph`i' = _se[exp0`shock'];
  gen segrech`i' = _se[rec0`shock']; 

  
  // Contracts
  
  ivreg2 F`i'.c $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bclinh`i' = _b[$`shock'linshock];  
  gen seclinh`i' = _se[$`shock'linshock]; 

  ivreg2 F`i'.c exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen bcexph`i' = _b[exp0`shock']; 
  gen bcrech`i' = _b[rec0`shock'];
  
  gen secexph`i' = _se[exp0`shock'];
  gen secrech`i' = _se[rec0`shock']; 
  
  
  ivreg2 F`i'.w $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bwlinh`i' = _b[$`shock'linshock];  
  gen sewlinh`i' = _se[$`shock'linshock]; 

  ivreg2 F`i'.w exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen bwexph`i' = _b[exp0`shock']; 
  gen bwrech`i' = _b[rec0`shock'];
  
  gen sewexph`i' = _se[exp0`shock'];
  gen sewrech`i' = _se[rec0`shock']; 
  
    
  ivreg2 F`i'.n $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bnlinh`i' = _b[$`shock'linshock];  
  gen senlinh`i' = _se[$`shock'linshock]; 

  ivreg2 F`i'.n exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen bnexph`i' = _b[exp0`shock']; 
  gen bnrech`i' = _b[rec0`shock'];
  
  gen senexph`i' = _se[exp0`shock'];
  gen senrech`i' = _se[rec0`shock']; 
  
};

  replace sumliny = bylinh`i' + sumliny;
  replace sumling = bglinh`i' + sumling;
  replace sumlinc = bclinh`i' + sumlinc; 
  replace sumlinw = bwlinh`i' + sumlinw; 
  replace sumlinn = bnlinh`i' + sumlinn; 
  
  replace sumexpy = byexph`i' + sumexpy;
  replace sumexpg = bgexph`i' + sumexpg;
  replace sumexpc = bcexph`i' + sumexpc; 
  replace sumexpw = bwexph`i' + sumexpw; 
  replace sumexpn = bnexph`i' + sumexpn; 
  
  replace sumrecy = byrech`i' + sumrecy;
  replace sumrecg = bgrech`i' + sumrecg;
  replace sumrecc = bcrech`i' + sumrecc; 
  replace sumrecw = bwrech`i' + sumrecw; 
  replace sumrecn = bnrech`i' + sumrecn; 
  
   gen multlinh`i' = sumliny/sumling;
   gen multexph`i' = sumexpy/sumexpg;
   gen multrech`i' = sumrecy/sumrecg;
  
  foreach var in bylin byexp byrec bglin bgexp bgrec bclin bcexp bcrec bwlin bwexp bwrec bnlin bnexp bnrec multlin multexp multrec {;
  
    quietly replace `var' = `var'h`i' if h==`i';
	
  };
  
  foreach var in ylin glin clin wlin nlin yexp gexp cexp wexp nexp yrec grec crec wrec nrec {;
  
    quietly replace up95b`var' = b`var'h`i' + 1.96*se`var'h`i' if h==`i';
	quietly replace lo95b`var' = b`var'h`i' - 1.96*se`var'h`i' if h==`i';
	
  };

};



forvalues i = 1/20 {;


ivreg2 F`i'.y $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bylinh`i' = _b[$`shock'linshock];  
  gen seylinh`i' = _se[$`shock'linshock];
  
ivreg2 F`i'.g $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bglinh`i' = _b[$`shock'linshock];  
  gen seglinh`i' = _se[$`shock'linshock]; 
  
 ivreg2 F`i'.c $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bclinh`i' = _b[$`shock'linshock];  
  gen seclinh`i' = _se[$`shock'linshock]; 
  
 ivreg2 F`i'.w $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bwlinh`i' = _b[$`shock'linshock];  
  gen sewlinh`i' = _se[$`shock'linshock]; 

 ivreg2 F`i'.n $`shock'linshock $`shock'linxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);

  gen bnlinh`i' = _b[$`shock'linshock];  
  gen senlinh`i' = _se[$`shock'linshock]; 

ivreg2 F`i'.y exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen byexph`i' = _b[exp0`shock'];
  gen byrech`i' = _b[rec0`shock'];
  
  gen seyexph`i' = _se[exp0`shock'];
  gen seyrech`i' = _se[rec0`shock']; 

ivreg2 F`i'.g exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen bgexph`i' = _b[exp0`shock']; 
  gen bgrech`i' = _b[rec0`shock'];
  
  gen segexph`i' = _se[exp0`shock'];
  gen segrech`i' = _se[rec0`shock'];
  
ivreg2 F`i'.c exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen bcexph`i' = _b[exp0`shock']; 
  gen bcrech`i' = _b[rec0`shock'];
  
  gen secexph`i' = _se[exp0`shock'];
  gen secrech`i' = _se[rec0`shock'];
  
ivreg2 F`i'.w exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen bwexph`i' = _b[exp0`shock']; 
  gen bwrech`i' = _b[rec0`shock'];
  
  gen sewexph`i' = _se[exp0`shock'];
  gen sewrech`i' = _se[rec0`shock'];
  

ivreg2 F`i'.n exp0`shock' rec0`shock' $`shock'nlxlist if L`p'.`omit'==0 & F`i'.`omit'==0 & `omit'==0, robust bw(auto);	

  gen bnexph`i' = _b[exp0`shock']; 
  gen bnrech`i' = _b[rec0`shock'];
  
  gen senexph`i' = _se[exp0`shock'];
  gen senrech`i' = _se[rec0`shock'];
  
  
  replace sumliny = bylinh`i' + sumliny;
  replace sumling = bglinh`i' + sumling;
  replace sumlinc = bclinh`i' + sumlinc;
  replace sumlinw = bwlinh`i' + sumlinw; 
  replace sumlinn = bnlinh`i' + sumlinn; 
  
  
  replace sumexpy = byexph`i' + sumexpy;
  replace sumexpg = bgexph`i' + sumexpg;
  replace sumexpc = bcexph`i' + sumexpc;
  replace sumexpw = bwexph`i' + sumexpw; 
  replace sumexpn = bnexph`i' + sumexpn; 
  
  replace sumrecy = byrech`i' + sumrecy;
  replace sumrecg = bgrech`i' + sumrecg;
  replace sumrecc = bcrech`i' + sumrecc;
  replace sumrecw = bwrech`i' + sumrecw; 
  replace sumrecn = bnrech`i' + sumrecn; 
  
  gen multlinh`i' = sumliny/sumling;
  gen multexph`i' = sumexpy/sumexpg;
  gen multrech`i' = sumrecy/sumrecg;
  
  foreach var in bylin byexp byrec bglin bgexp bgrec bclin bcexp bcrec bwlin bwexp bwrec bnlin bnexp bnrec multlin multexp multrec seyexp seyrec segexp segrec secexp secrec {;
  
    quietly replace `var' = `var'h`i' if h==`i';
	
  };
  
  foreach var in ylin glin yexp gexp yrec grec clin cexp crec wlin wexp wrec nlin nexp nrec {;
  
    quietly replace up95b`var' = b`var'h`i' + 1.96*se`var'h`i' if h==`i';
	quietly replace lo95b`var' = b`var'h`i' - 1.96*se`var'h`i' if h==`i';
	
  };

  
};

display as text "MULTIPLIERS:  2 STEP";

rename multlin multlin2;
rename multexp multexp2;
rename multrec multrec2;

outsheet h multlin2 multexp2 multrec2 using junk2step.csv if h<=20, comma replace;

outsheet h byexp byrec bgexp bgrec bcexp bcrec seyexp seyrec segexp segrec using junk2stepirfs.csv if h<=20, comma replace;

label var bglin "Gov, linear model";
label var bylin "GDP, linear model";
label var bgexp "GOV, expansion";
label var byexp "GDP, expansion";
label var bgrec "Gov, recession";
label var byrec "GDP, recession";

label var bclin "Contracts, linear model"; 
label var bwlin "Wages, linear model"; 
label var bnlin "Non-Wage Missing, linear model"; 


tw (rarea up95bglin lo95bglin h, bcolor(gs12) clw(medthin medthin)) 
  (scatter bglin h, c(l ) clp(l ) ms(i ) clc(black) mc(black) clw(medthick)) if h<=20,
  saving(../output/ramey_zubairy_replication/junkg.gph,replace);

tw (rarea up95bylin lo95bylin h, bcolor(gs12) clw(medthin medthin)) 
  (scatter bylin h, c(l ) clp(l ) ms(i ) clc(black) mc(black) clw(medthick)) if h<=20,
  saving(../output/ramey_zubairy_replication/junky.gph,replace);
  
tw (rarea up95bclin lo95bclin h, bcolor(gs12) clw(medthin medthin)) 
  (scatter bclin h, c(l ) clp(l ) ms(i ) clc(black) mc(black) clw(medthick)) if h<=20,
  saving(../output/ramey_zubairy_replication/junkc.gph,replace);
  
tw (rarea up95bwlin lo95bwlin h, bcolor(gs12) clw(medthin medthin)) 
  (scatter bwlin h, c(l ) clp(l ) ms(i ) clc(black) mc(black) clw(medthick)) if h<=20,
  saving(../output/ramey_zubairy_replication/junkw.gph,replace);
  
tw (rarea up95bnlin lo95bnlin h, bcolor(gs12) clw(medthin medthin)) 
  (scatter bnlin h, c(l ) clp(l ) ms(i ) clc(black) mc(black) clw(medthick)) if h<=20,
  saving(../output/ramey_zubairy_replication/junkn.gph,replace);
  
tw (rarea up95bgrec lo95bgrec h, bcolor(gs12) clw(medthin medthin))
    (scatter up95bgexp bgexp lo95bgexp bgrec h, clw(medthin medthick medthin medthick)
  c(l l l l l) clp(- l - l) clc(red red red blue ) ms(i o i i i i) mc(red red red blue)) if h<=20,
  saving(../output/ramey_zubairy_replication/junkgnl.gph,replace);

tw (rarea up95byrec lo95byrec h, bcolor(gs12) clw(medthin medthin))
    (scatter up95byexp byexp lo95byexp byrec h,  clw(medthin medthick medthin medthick)
  c(l l l l l) clp(- l - l) clc(red red red blue ) ms(i o i i i i) mc(red red red blue)) if h<=20,
  saving(../output/ramey_zubairy_replication/junkynl.gph,replace);
  
  
tw (rarea up95bcrec lo95bcrec h, bcolor(gs12) clw(medthin medthin))
    (scatter up95bcexp bcexp lo95bcexp bcrec h,  clw(medthin medthick medthin medthick)
  c(l l l l l) clp(- l - l) clc(red red red blue ) ms(i o i i i i) mc(red red red blue)) if h<=20,
  saving(../output/ramey_zubairy_replication/junkcnl.gph,replace);


outsheet h *lin up95* lo95* if h <= 20 using ../output/ramey_zubairy_replication/plotvars_BP.csv, comma replace;
  
   
graph combine ../output/ramey_zubairy_replication/junkg.gph ../output/ramey_zubairy_replication/junkc.gph ../output/ramey_zubairy_replication/junkw.gph ../output/ramey_zubairy_replication/junkn.gph, cols(2) iscale(0.5);

graph export "../output/ramey_zubairy_replication/combo_`shock'.pdf", replace;



capture log close;

 
